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Abstract. We present here in full details the linear secular theory with tidal damping that was used to constraint the fit of the 
HD10180 planetary system in ( [Lovis et al.|201 l| l. The theory is very general and can provide some intuitive understanding of 
the final state of a planetary system when one or more planets are close to their central star. We globally recover the results of 
(|Mardling'2007'l, but we show that in the HD209458 planetary system, the consideration of the tides raised by the central star 
on the planet lead to believe that the eccentricity of HD209458b is most probably much smaller than 0.01. 



1. Introduction 

In several planetary systems, some planets are very close to 
their central star and thus subject to strong tidal interaction. If 
a planet were alone around its star, this will led to a circulariza- 
tion of its orbit. However, for multi-planet system, due to the 
secular interaction between the planets, the final evolution of 
the system may be different, with residual eccentricity (Wu &' 



Goldreich,2002 , Mardling 2007t|Batygin et al.|2009 ; Mardhng 



2010 Lovis et al.|201 1 1. This is also why fitting a circular orbit 
to the innermost planets in a system subject to tidal dissipation 
will not insure that its eccentricity remains small, as the secular 
interactions may drive it to large values (Lovis et al. 201 1\ . One 
thus needs to take into account the fact that the observed sys- 
tem is the result of a tidal process ( Lovis et al.|201 1 1. Here, we 
develop in full details the method that has been used in ( |Lovis| 
|et al.|20lT I for the fit to a tidally evolved system. The theory is 



very general and is compared to previous results of ( [Mardhng 
[20071 [20T0l l. 



tively the eccentricity and the longitude of the periastron of the 
k-th planet, the secular equations read 



dt 



[z] = iA [z] , with [z] 



Zn ) 



(1) 



and where A is a real matiix whose elements are ( |Laskar &[ 
|Robutel|1995| ) 



fef '«o \ajj jf-^, moflt \akj 



2„^.!!?i^C,p if j<k, 
moat \akj 



(2) 



liij — C2 — 
mo \aj 



if j>k. 



2. Model 

2.1. Newtonian interaction 

Without mean-motion resonances, the long term evolution of 
a conservative multiplanetary system is given, in first order. 



by the Laplace-Lagrange linear secular equations (see Laskar 



[1990 ). With this approximation, inclinations and eccentricities 
are decoupled and follow the same kind of evolution. For sim- 
phcity, in this letter we will focus only on coplanar systems. 
Let n be the number of planets. Using the classical complex 
variables Zk = e^e"^' , k - 1, . . . , n, where and mk are respec- 



In these expressions, nik and Ok, are the mass and semi-major 
axis of the A:-th planet, while the mean motion is defined 
by n^fl^ = G(mo + nik). By assumption, planets are ordered by 
increasing semi-major axis while the index stands for the star. 
The functions C2(a) and C^ia) are defined by mean ofbf ,the 
usual Laplace coefficients (e.g. Laskar & Robutel|1995 1, as: 



3*"( 

-'3/2'- 



(3) 



1 



C3(a) = jah\)>^(a) . 
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2.2. Other effects 

The above secular equations ([T]) describe only the Newtonian 
interactions between point mass planets. In order to study ex- 
oplanetary systems with short period planets, it is often nec- 
essary to add corrections due to relativity, oblateness and tidal 
friction, at least to the inner planets. The spatial secular equa- 
tions of motion resulting from all these effects are given in 



(|Lambeck"1980',' Eggleton & Kiseleva-Eggleton|200T] [Ferraz- 
^ello et al. 2008). To first order in eccentricity, and in the pla- 
nar case, they modify the diagonal terms of the matrix A in 
two different ways. There are conservative terms that are purely 
imaginary, and dissipative ones which are real (and negative). 
We thus define two new diagonal matrices 6A - 2/=i,5 '5^*'' 
6B = 2 1=4,5 such that the full secular evolution is given by 



:T-M = (iAtot-5B)k], 

at 



(4) 



with Atot - A + 6A. The effect of relativity on the A:-th planet is 
conservative, and in first order in eccentricity, it leads to 



(5) 



The effect of the oblateness of bodies generated by their proper 
rotation are 



and 

,(3) 
kk 



k2,k nk I mo + Mk \lRk\^ l<^k^^ 



Uk \nk 



6A 



k2,k rik I mo + Mk \ / \^ / Wo ^ ^ 



mo 



(6) 



(7) 



for the oblateness of the A:-th planet and the oblateness of the 
star. A:2_t, w^, and R}. are respectively the second Love number, 
the proper rotation rate, and the radius of the ^-th body. Tidal 
effects have two contributions. With the same approximation, 
we have 



11 

18 nj a ' 



with 



mo \ / Rk 



Kk = k2,knk 

■ mi l\ak, 

for the tides raised on the A:-th planet by the star, and 

1 1 Wo \ K[ 



^41' =2711 



18 rikj Qo 



with 



mu \(Ro 



(8) 



(9) 



(10) 



(11) 



for the tides raised on the star by the A;-th planet. We consider 
here the "viscous" approach (Singer 1968 Mignard 1979| l, 
where the quality factor of the A:-th body is Qk = {nk{ts.t)ky^ 
and {At)k is a constant time lag. 



3. Resolution 

3.1. Conservative case 

When there is no dissipation (6B - 0) the system is classi- 
cally resolved by diagonalizing the matrix Atot through a linear 
transformation 

[z] = SoM . (12) 
In the new variables, the equations of motion become 

d _i 

— [m] = iDo[u] , where Do = Atot^o (13) 

is the diagonal matrix diag(gi, . . . ,g„) of the eigenvalues gk of 
Atot. We have then 



Uk(t) = M;t(0)e' 



igtt 



(14) 



Each proper mode Uk describes a circle in the complex plane at 
a constant frequency and with the radius |Mt(0)|. The evolu- 



tion of the planetary eccentricities are then given by ( 12 1. They 
are linear combinations of the proper modes. 

The only differences between A and Atot are in the diago- 
nal terms. Those of Atot are larger or equal to those of A. As 
a consequence, using Atot instead of A makes the fundamen- 
tal frequencies g^ larger and the coupling between the proper 
modes lower (the evolution of the eccentricity of each planet 
is almost given by one single proper mode, the other modes 
generate only small oscillations). 

3.2. General solution 

In the full linear secular equation Q, the matrix that has to 
be diagonalized is now iAtot - SB, where the dissipation part 
6B comes only from tides. In general, the elements of 6B are 
much smaller than those of the diagonal of Atot. SB will thus 
be considered as a perturbation of the conservative evolution 
given by Atot. Let 



S =5o(l +i(55i) 



(15) 



be the matrix of the linear transformation that diagonalizes the 
full system. As 6B is a perturbation of Atot, we will make the 
hypothesis that the matrix 6S i is also a perturbation of the ma- 
trix ^o. At first order, the inverse of S is 



S-^ = (1 -i55i)5o' 



and the new diagonal matrix is D = iDo - 6Di, with 

6D, ^S^,\6B)So-[6SuDo\ , 

where the bracket is defined by [6S i,Do] - 6S ]Do 
For 6Di to be actually diagonal, 6S i is given by 



(SSi)jk 



1 



-{So\6B)So) 



jk 



(16) 

(17) 
DoSSi. 

(18) 



8k - 8j 

As Do is diagonal, all terms in the diagonal of [6S i , Dq] 
vanish. Thus, the diagonal terms of 6S i do not appear in the 



computation of 6Di ( 17 1 and they can be set equal to zero. Let 
6Di = diag(yi, . . . ,j„). From ( [T7| , we have then 



7k 



^{s,U6B)Sol. 



(19) 
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These jk are real and positive. It turns out that the imaginary Table 1. Data for HD 209458b. 



part of D is still the one of the conservative case Do (13 i. The 
proper frequencies git are not affected by the dissipation 6B. 
However, each proper mode now contains a damping factor 
given by ( [T9] l. The equations of motion in the new variables 
now read 



dt 



[u] = diagfei - n, . . . , ig„ - yn)[u] 



(20) 



and the solutions are 
Uk(t) = Mi(0)e->'"e'«" , 



It should be stressed that even if only one planet undergoes 
tidal dissipation (only {6B)[i is different from for example), 
because of the linear transformation ^o, all the eigenmodes can 



be damped ( 19 1 



4. Two planet case 

In a simpler two planet system where only the first one under- 
goes tidal friction, 6B = diag(y, 0), the two proper frequencies 
are given by 



82 



^ (r + Vr2 - 4a) 
^ (r - Vr2 - 4a) 



(22) 



where T and A are the trace and determinant of A,„,. From (19\ , 
it can be shown that the two dissipation factors are 



ri 



72 



1 + 



122 



1 



^1 -82 

An -A22 
81 -82 



7 ■ 



7 ■ 



(23) 



The sum yi + 72 is equal to y. There is thus always one eigen- 
mode damped in a timescale shorter than 2y ' while the other 
is damped in a timescale larger than 2y ' . In the particular case 
where An = A22, we have yi = 72 = y/2. 

Once the first eigenmode is damped, the ratio between the 
two eccentricities and the difference between the two longi- 



tudes of periastron are deduced from ( 12 1. We have 




tzT] - — , if yi > y2 



m\ — tiT2 — TT , if yi < y2 



(24) 



It should be noted that the matrix 6S 1 introduces small cor- 
rections in the difference between the longitudes of periastron 
which are not taken into account in (|24]). 

5. Application to HD 209458b 

Here we compare the results of this paper with those of 
Mardhng|([2007ll on the example of HD 209458b (Table[T). As 
m ( |Mardling||2007l ), we first assume that the non zero eccen- 
tricity of this planet is due to the presence of a m2 - 0.1 Mj 



HD 209458b 



Period (day) 


3.5247 


mo (Mo) 


1.10 


nil (Mj) 


0.64 ± 0.06 


a I (AU) 


0.045 




0.014 ±0.009 


Ri (Rj) 


1.32 ±0.03 



note: As in l |Mardling|2007[ l, all parameters come from ( Burrows et al.| 
(21) |2007| l except the eccentricity that comes from ( ,Laughlin et al.,2005| l. 



companion at 02 - 0.4 AU with an eccentricity 62 = 0.4. For 
this study, eccentricities are large and modify the frequencies 
given by the analytical expression of the matiix (|4]i. Thus, 
we chose to compute the matrix Atot using a frequency analy- 
sis on a numerical integration of the system without dissipation 
exact in eccentricity and expanded up to the 4 order in the ratio 
of the semi-major axes (e.g. Mardling & Lin|[2002 Laskar & 



Boue )2()T0| ). At first order, the eccentricity variables zi and Z2 
are linear combinations of two eigenmodes mi and U2 (21 1 



Zl(f) ^SuUiit) + S 12M2(0 , 
Z2(t) = 5 21 Ml (f) + S22U2(t) , 



where S is given by (15 1. With Qi - 

-1 



(25) 



10^ and wi = mi, the 
- y-i =46 Myr and 



two damping timescales (23 1 are 
yj ^ = 589 Gyr. With an age estimate of 5.5 Gyr for this system 
([Burrows et al.|[2007 1, the first eigenmode should be damped 
and the modulus of the second should remain almost constant. 
In consequence, both eccentricity variables should be propor- 
tional to U2. Their modulus should thus be constant and verify 
([24)1, or equivalently, ei ~ 62^5 12/^22 = 0.0025. In (Mardlmgl 
2007 Fig. 3), this value is larger, ei 



0.012. The difference 
comes from the tidal deformation of the planet that leads to the 
coefficient dAf^ in Eq. ([sj). This was not taken into account 



in (Mardling 20071, and it accelerates the precession of the 
periapse of HD 209458b by a factor 6.6 (see Fig. [T^). In fig- 
ure [T^, the initial Q-\alne of the planet is set artificialy to 100 



to enable a direct comparison with the figure 3 of (Mardling 
|2007[ ). As said by Mard ling] ( |2007| l, and showed in this paper, 
the 2-value only affects the damping timescales but not the pre- 
cession frequencies, nor the eccentricity amplitudes. However, 
with a larger precession frequency, the matrix Atot is closer to 
a diagonal matrix. The two planets are less coupled and the ra- 
tio S12/S22 ( |25] l, equal to the final eccentricity ratio ei/e2, is 
smaller. 

One way to recover the final eccentiicity of HD 209458b 
is to increase the mass of the companion up to m2 - 0.608My 
(Fig. [TJ)). Here, our aim is not to explain the large eccentric- 
ity of HD 209458b, but simply to illustrate the results of the 
section [3?2l 

As the precession of the periastron of the inner planet is 
faster than in (|Mardling 20071 Fig. 3), we decreased the ini- 



tial Q-value to 15.15 to accelerate the damping and to obtain 
an evolution with the same gi/j\ ratio as in ( |Mardling||2007| l 
(Fig.[T];). As said before, this does not change the final eccen- 
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10 20 30 40 50 
time (10* yr) 




1 = 15.15 
m2 = 0.608Mj 



1 2 3 4 5 6 7 
time (10* yr) 




0.12 
0.1 
0.08 
0.06 
0.04 
0.02 



10 20 30 40 
time (10* yr) 



50 



(d) 

Qi = 15.15 
m2 = 0.608Mj 
with a planet 




10 20 30 40 
time (10* yr) 



50 



Fig. 1. Tidal effects on the eccentricity of HD 209458b with one or 2 companions, (a) The companion is a 0. IMy planet at 0.4 AU 
with 62 - 0.4 as in ( Mardling|2007 1. The blue curve has been obtained without considering the conservative effect of the tides 
(Mf/in eq.js}. Thered curve in (a), and all the evolutions in the subfigures (b), (c), and (d) take into account this effect, (b) The 
mass of the companion is set to 0.608My in order to recover the final excentricity of Mardling's simulation. The red curve is the 
result of a numerical integration of the full secular equations exact in eccentricity. The green one is the analytical solution of the 
linearized problem, (c) Same as (b) except for the initial Q-yalne of HD 209458b which is set to 15.15 to enable the visualisation 
of both the damping and the oscillation of the eccentricity, (d) Same as (c) with an additional O.lMy companion at 1.0 AU with 
63 =0.1. 



tricity, but it illustrates better the damping of the first mode with 
frequency gi =0.14 deg/yr. 

After the damping of the first eigenmode, the eccentrici- 
ties are not oscillating because there remains only one eigen- 
mode with a non-zero amplitude. Both eccentricity variables z\ 
and Z2 describe a circle in the complex plane at the same fre- 
quency g2- But if a third planet is added to the system, a new 
eigenmode appears with a frequency g^- Then eccentricities are 
oscillating (Fig. [T|J). It should be noted that a relative inclina- 
tion between planets can also generate an other eigenmode and 
make eccentricities oscillate fMardling'2010'). However, in the 
linear approach eccentricities and inclinations are decoupled. It 
is thus necessary to have large eccentricities or inclinations to 
have significant oscillations. 

We now wonder which companion parameters can lead to 
an eccentricity ei - 0.01 for HD 209458b. As the system con- 
tains two planets, eccentricities are at most combination of two 
eigenmodes. But since y"' - 46Myr is less than the age of 
the system (5.5 Gyr), at least one of the eigenmode is damped. 



However both eigenmodes cannot have zero amplitude, else the 
two orbits would be circular. Thus, let us assume that remains 
only a single eigenmode, the one with the longer damping 
timescale. Then, given a semi-major axis 02 and a mass 1112, the 



eccentricity of the companion is obtained through ( 24 1 at first 
order. In practice we integrated numerically the system without 
dissipation, and found the eccentricity 62 that cancels the am- 
plitude of the rapidly damped eigenmode. Results are shown 
with solid curves in figure |2^ and|2|5. Once the current eccen- 
tricity 62 is given, the initial value (5.5 Gyr ago) is estimated 
assuming an exponential decay with a damping factor given by 
( [23] l (see the dotted curves Fig. [2^ and Fig.|2]5). The frequencies 
gi- and the coefficients A^^ were obtained numerically using a 
frequency analysis. Parameters leading to initial eccentricities 
larger than 1 are excluded. They correspond to the grey regions 
in figure[2] Although planets are less coupled than in (Mardling 
|2007[ ), there is still a large range of initial conditions leading 
to a state compatible with = 0.01. However, the stellar re- 
flex velocity due to the companion at periastron (Fig. |2}; and 
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Fig. 2. Eccentricity €2 of the hypothetical companion of HD 209458b with e\ - 0.01 assuming that the eigenmode with the 
shortest dissipation timescale is damped (black curves in panel a) and b)). a) The mass of the companion is fixed to OT2 - Q.2Mj. 
Negative values of 62 correspond to Am - 180 deg while positive ones mean Am - Odeg. The dotted line is the eccentricity that 



the companion would have had 5.5 Gyr ago assuming a dissipation factor computed with (23 1. b) Same as a) for different masses 
m2 while the semi-major axis is fixed and set to 02 - 0.4 AU. c) Stellar reflex velocity due to the companion at periastron with 
the eccentricity of the figure a), d) idem for the eccentricity of the figure b). In grey regions, the eccentricity of the companion 
should have been larger than 1 in the past. The configuration appearing in all panels with the same orbital parameters is marked 
by a fill circle. 



[2]l) is above the detectability threshold of about 3 m.s For 
example, with 02 = 0.25 AU, and m2 - Q.QSMj, the current 
eccentricity is 62 = 0.34 and the maximal stellar reflex velocity 
Vq = 3.9 m.s"'. It thus seems that such a planet cannot exist. 
Indeed, observations do not constrain strongly the eccentricity 



of HD 209458b and a circular orbit is not ruled out (Laughlin 
let al.|2005| l. 

6. Orbit fitting : the I-ID10180 case 

The analysis of the radial velocities measures of HD10180 re- 
vealed the potential existence of 7 planets in this system (Lovis' 



|et al.|,2011j . The innermost planet, HD10180b, is a terrestrial 
planet (nihsini - 1.35Me) with period of x 1.177 days and 
semi-major axis ai, = 0.0223 AU. The planet is thus subject to 
strong tidal interaction with the central star. During the first fit 
( |Lovis et al.|201 1| Table 3) , it was thus assumed that its eccen- 
tricity has been damped to very small values, and its value was 
fixed to et = 0. Nevertheless, if the system is then numerically 
integrated over 250 kyr (Fig. |4] (red curve)), due to secular in- 
teractions with the other planets, e^, grows very rapidly to high 
values, reaching nearly 0.9. 




0.4 0.6 
time (Gyr) 

Fig. 3. Tidal evolution of the amplitude of the proper modes 
|mi| (red), \u2\ (green), \ui\ (blue), and U4 (pink) resulting from 
the tidal dissipation on planet HDIOI8OZ7 with fca/fi - 0.0015 
( |Lovis et aLpOTTj i. 



When general relativity (GR) is included in the numerical 
integration, the main effect is to increase the diagonal terms of 
the secular matrix (Eq|5]). As a result, the secular variations of 
eh are much smaUer (Fig.|4](green curve)), but still reach 0.2. 
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and the additional contribution ( 26 1 can then be computed in 



100 150 
time (kyr) 



250 



the fitting process. Practically, in an iterative fit taking into 
account the Newtonian interactions, the transformation matrix 
just needs to be computed once, or twice if one wants to 
recompute the matrix when the convergence to a final so- 
lution is obtained. In ( |Lovis et al.|201 l| l, the final values were 
Ml = 0.0017, M2 = 0.044 for R = 350, with a final ^Jx^ = 1.24, 
very close to the residuals obtained in absence of constraint ( 
V?=1.22). 

In this constrained solution, the initial value of e/, is still 
0, but the secular change due to planetary interactions is much 
smaller (Fig|4](blue curve)), which ensure a more stable behav- 
ior to the system. 



Fig. 4. Evolution of the eccentricity of planet HDlOlSOb over 
250 kyr starting with e^, = at f = (present time) for three 
different models : In red, the numerical integration is purely 
Newtonian and do not take into account general relativity (GR). 
In green, GR is taken into account in the integration. In blue, 
GR is taken into account and the fit is made with the tidal dis- 



sipation constraint (26 1. 



The strategy that was adopted for the final fit of ( Lovis et al. 
201 1) was to include in the fit the constraint that the planetary 
system that is searched for is the result of the tidal evolution, as 



described in section 3.2 As the planet has a mass comparable 
to the Earth, it can be assumed to be terrestrial, and thus to have 
a dissipation factor of the same order of magnitude as (or larger 
than) Mars ^2/6 = 0.0015, which is the smallest value among 
the terrestrial planets in the Solar System. The damping fac- 



tors e"''' t can thus be computed through (Eq. 19 1 for all proper 
modes Uf, ( [Lovis et al.|20ri^ Table 5). The resulting dampings 
of the amplitudes of the proper modes Uk are given in Fig|3] 

From this computation, as the age of the system is estimated 
to be of about 4 Gyr, it can be seen that the first two proper 
modes amplitudes mi and U2 are certainly reduced to very small 
values. If the damping factor k2/Q were 10 times smaller, the 
conclusion would be nearly the same, as the only change in 
Figj3] would be a change in the time scale of the figure, the 
units being now 10 Gyr instead of 1 Gyr. 

In order to include the constraint on the tidal damping in 
the fit, one can then add to the minimization the additional 
term 



-.R(\ui\^ + \u2f) 



(26) 



where R is an empirical constant that needs to be set to a value 
that will equilibrate the additional constraint with respect to the 
value of the in absence of constraint. After various trials, 
R = 350 was used in ( |Lovis et al.|201 l| l. 

The computation of the ampUtude of the proper modes 
\uk\ during the fit is made iteratively. Once a first orbital so- 
lution is obtained, the Laplace-Lagrange linear system (Eq{T]i 
is computed and thus also the matrix So of transformation to 
proper modes (EqfT2|. For a given initial condition (zk) ob- 
tained through the fit, the proper modes are computed with 



(27) 



7. Conclusion 

We have presented here in full details the secular theory with 
tidal dissipation that was outlined in ( Lovis et al.|201 1 1 for the 
system HD 10180. The use of Lagrange-Laplace linear theory 
can include very easily the linear contribution due to tidal dissi- 
pation and provide an intuitive background for studying multi- 
planetary systems when one or several planets are close to their 
central star and subject to tidal damping. Although we have 
hmited here the study to the planar case, this formalism can be 
easily extended to mutually inclined systems. 

For the system HD209458, we could retrieve globally the 
results of ( MardUng|2007 '), although we find that a companion 
with mass m2 = 0.1 My with a2 = 0.4 AU and 62 = 0.4 will not 
lead to ei = 0.012, but to a much smaller value of ei = 0.0025. 
This is due to the additional tides raised by the star on the planet 
(Eqjij) in the contribution to the secular equations (Eq|4]). 

We have examined other configuration that could lead to a 
final eccentricity ei > 0.01 for HD209458b, but our conclu- 
sion is negative, as we found that a potential companion, large 
enough to lead to a final eccentricity ei > 0.01 leads to suf- 
ficiently large stellar motion that it should have akeady been 
detected, assuming a detectability threshold of 3 m.s Our 
conclusion is thus that the most probable outcome is that the 
actual eccentricity of HD209458b has a much smaller value 
than 0.01. 
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